home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / libgutil / geom.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  11KB  |  496 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. /*
  18.  *    geom -
  19.  *        A software geometry package for computer graphics.
  20.  *
  21.  *                David J. Brown - 1983 
  22.  */
  23. #include "values.h"
  24. #include "math.h"
  25. #include "stdio.h"
  26.  
  27. #define STACKSIZE    20
  28.  
  29. static float Vm_Stack[STACKSIZE][4][4];
  30. static int   SP = -1;
  31. static float Vm[4][4];    /* Viewing matrix */
  32. static float Vsx, Vsy, Vcx, Vcy;    /* Viewport parameters */
  33.  
  34. /*
  35.  *     initmatstack - 
  36.  *        Initialize the matrix stack and the viewport
  37.  */
  38. i_initmatstack()
  39. {
  40.     i_initmatrix();
  41.     i_initviewport();
  42. }
  43.  
  44. /*
  45.  *     initviewport - 
  46.  *        Initialize the viewport
  47.  */
  48. i_initviewport()
  49. {
  50.     Vsx = 1.0;
  51.     Vsy = 1.0;
  52.     Vcx = 0.0;
  53.     Vcy = 0.0;
  54. }
  55.  
  56. /*
  57.  *     initmatrix - 
  58.  *        Set the matrix to the identity    
  59.  */
  60. i_initmatrix()
  61. {
  62.     i_initmat(Vm);
  63. }
  64.  
  65. i_initmat(matrix)
  66. float *matrix;
  67. {
  68.     *matrix++ = 1.0;    /* row 1    */
  69.     *matrix++ = 0.0;
  70.     *matrix++ = 0.0;
  71.     *matrix++ = 0.0;
  72.     *matrix++ = 0.0;    /* row 2    */
  73.     *matrix++ = 1.0;
  74.     *matrix++ = 0.0;
  75.     *matrix++ = 0.0;
  76.     *matrix++ = 0.0;    /* row 3    */
  77.     *matrix++ = 0.0;
  78.     *matrix++ = 1.0;
  79.     *matrix++ = 0.0;
  80.     *matrix++ = 0.0;    /* row 4    */
  81.     *matrix++ = 0.0;
  82.     *matrix++ = 0.0;
  83.     *matrix++ = 1.0;
  84. }
  85.  
  86. /*
  87.  *      viewport - 
  88.  *        Set the current viewport
  89.  */
  90. i_viewport( l, r, b, t)
  91. float l, r, b, t;
  92. {
  93.     Vsx = r-l;
  94.     Vsy = t-b;
  95.     Vcx = (r+l)/2.0;
  96.     Vcy = (t+b)/2.0;
  97. }
  98.  
  99. /*
  100.  *      pushmatrix - 
  101.  *        Push the current matrix onto the stack
  102.  */
  103. i_pushmatrix()
  104. {
  105.     int row, col;
  106.  
  107.     if (SP == STACKSIZE-1)
  108.         puts("Stack is Full");
  109.     else {
  110.     SP++;
  111.         for(row=0 ; row<4 ; row++)
  112.             for(col=0 ; col<4 ; col++)
  113.                 Vm_Stack[SP][row][col] = Vm[row][col];
  114.     }
  115. }
  116.  
  117.  
  118. /*
  119.  *      popmatrix - 
  120.  *        Pop the current matrix off the stack
  121.  */
  122. i_popmatrix()
  123. {
  124.     int row, col;
  125.  
  126.     if (SP == -1)
  127.         puts("Stack is Empty");
  128.     else {
  129.         for(row=0 ; row<4 ; row++)
  130.             for(col=0 ; col<4 ; col++)
  131.                 Vm[row][col] = Vm_Stack[SP][row][col];
  132.     SP--;
  133.     }
  134. }
  135.  
  136. /*
  137.  *    scale - 
  138.  *        Scale the current matrix
  139.  */
  140. i_scale(Sx, Sy, Sz)
  141. float Sx, Sy, Sz;
  142. {
  143.     int i;
  144.  
  145.     for(i=0 ; i<4 ; i++)
  146.     Vm[0][i] *= Sx;
  147.     for(i=0 ; i<4 ; i++)
  148.     Vm[1][i] *= Sy;
  149.     for(i=0 ; i<4 ; i++)
  150.     Vm[2][i] *= Sz;
  151. }
  152.  
  153. /*
  154.  *    loadmatrix - 
  155.  *        Load the current matrix
  156.  */
  157. i_loadmatrix( m )
  158. float m[4][4];
  159. {
  160.     int i, j;
  161.  
  162.     for(i=0; i<4; i++) 
  163.     for(j=0; j<4; j++) 
  164.         Vm[i][j] = m[i][j];
  165. }
  166.  
  167. /*
  168.  *    translate - 
  169.  *        Translate the current matrix
  170.  */
  171. i_translate(Tx, Ty, Tz)
  172. float Tx, Ty, Tz;
  173. {
  174.     Vm[3][0] += (Tx*Vm[0][0] + Ty*Vm[1][0] + Tz*Vm[2][0]);
  175.     Vm[3][1] += (Tx*Vm[0][1] + Ty*Vm[1][1] + Tz*Vm[2][1]);
  176.     Vm[3][2] += (Tx*Vm[0][2] + Ty*Vm[1][2] + Tz*Vm[2][2]);
  177. }
  178.  
  179. /*
  180.  *    rotate - 
  181.  *        Rotate about an axis
  182.  */
  183. i_rotate(angle, axis)
  184. float angle;
  185. char axis;
  186. {
  187.     int row,col;
  188.     float temp[4];
  189.     float cosine, sine;
  190.  
  191.     for(col=0; col<4 ; col++)    /* init temp to zero matrix */
  192.         temp[col] = 0;
  193.  
  194.     angle = angle*(3.1415926535/180.0);
  195.     cosine = cos(angle);
  196.     sine = sin(angle);
  197.     switch(axis){
  198.     case 'x':    
  199.     case 'X':    
  200.         for(col=0 ; col<4 ; col++)
  201.             temp[col] = cosine*Vm[1][col] + sine*Vm[2][col];
  202.         for(col=0 ; col<4 ; col++) {
  203.         Vm[2][col] = - sine*Vm[1][col] + cosine*Vm[2][col];
  204.             Vm[1][col] = temp[col];
  205.     }
  206.         break;
  207.  
  208.     case 'y':
  209.     case 'Y':
  210.         for(col=0 ; col<4 ; col++)
  211.             temp[col] = cosine*Vm[0][col] - sine*Vm[2][col];
  212.         for(col=0 ; col<4 ; col++) {
  213.             Vm[2][col] = sine*Vm[0][col] + cosine*Vm[2][col];
  214.             Vm[0][col] = temp[col];
  215.         }
  216.     break;
  217.  
  218.     case 'z':
  219.     case 'Z':
  220.         for(col=0 ; col<4 ; col++)
  221.             temp[col] = cosine*Vm[0][col] + sine*Vm[1][col];
  222.         for(col=0 ; col<4 ; col++) {
  223.             Vm[1][col] = - sine*Vm[0][col] + cosine*Vm[1][col];
  224.             Vm[0][col] = temp[col];
  225.         }
  226.     break;
  227.     }
  228. }
  229.  
  230. /* 
  231.  *        multmatrix - 
  232.  *        Premultiply the current matrix 
  233.  */
  234. i_multmatrix(icand)
  235. float icand[4][4];
  236. {
  237.     int row, col;
  238.     float temp[4][4];
  239.  
  240.     for(row=0 ; row<4 ; row++) 
  241.         for(col=0 ; col<4 ; col++)
  242.             temp[row][col] = icand[row][0] * Vm[0][col]
  243.                            + icand[row][1] * Vm[1][col]
  244.                            + icand[row][2] * Vm[2][col]
  245.                            + icand[row][3] * Vm[3][col];
  246.     for(row=0 ; row<16 ; row++)
  247.          *(Vm[0]+row) = *(temp[0]+row);
  248. }
  249.  
  250. /*
  251.  *     transform -
  252.  *        Transform the coordinates using the current matrix.
  253.  */
  254. i_transform(x,y,z,tx,ty,tz)    
  255. float x,y,z;
  256. float *tx,*ty,*tz;
  257. {
  258.     int col;
  259.     float Vc[4];
  260.  
  261.     for(col=0 ; col<3 ; col++)
  262.         Vc[col] = x*Vm[0][col] + y*Vm[1][col] + z*Vm[2][col] + Vm[3][col];
  263.     *tx = Vsx*Vc[0]+Vcx;
  264.     *ty = Vsy*Vc[1]+Vcy;
  265.     *tz = Vc[2];
  266. }
  267.  
  268. /*
  269.  *     transform4 -
  270.  *        Transform the coordinates using the current matrix.
  271.  */
  272. i_transform4(x,y,z,w,tx,ty,tz,tw)    
  273. float x,y,z,w;
  274. float *tx,*ty,*tz,*tw;
  275. {
  276.     int col;
  277.     float Vc[4];
  278.  
  279.     for(col=0 ; col<4 ; col++)
  280.         Vc[col] = x*Vm[0][col] + y*Vm[1][col] + z*Vm[2][col] + w*Vm[3][col];
  281.     *tx = Vsx*Vc[0]+Vcx;
  282.     *ty = Vsy*Vc[1]+Vcy;
  283.     *tz = Vc[2];
  284.     *tw = Vc[3];
  285. }
  286.  
  287.  
  288. /*
  289.  *     polarview -
  290.  *        Concatenate a polarview matrix 
  291.  *
  292.  */
  293. i_polarview(dist, azimuth, incidence, twist)
  294. float dist;
  295. float azimuth, incidence, twist;
  296. {
  297.     i_translate(0.0, 0.0, -dist);
  298.     i_rotate(-twist,'z');    
  299.     i_rotate(-incidence,'x');
  300.     i_rotate(-azimuth,'z');
  301. }
  302.  
  303. /*
  304.  *     lookat -
  305.  *        Concatenate a lookat matrix 
  306.  */
  307. i_lookat(vx, vy, vz, px, py, pz, twist)
  308. float vx, vy, vz, px, py, pz;
  309. float twist;
  310. {
  311.     float sine, cosine, hyp, hyp1, dx, dy, dz;
  312.     float mat[4][4];
  313.  
  314.     i_rotate(-twist,'z');
  315.     dx = px - vx;
  316.     dy = py - vy;
  317.     dz = pz - vz;
  318.     hyp = dx * dx + dz * dz;    /* hyp squared    */
  319.     hyp1 = sqrt(dy*dy + hyp);
  320.     hyp = sqrt(hyp);        /* the real hyp    */
  321.     i_initmat(mat);
  322.     if (hyp1 != 0.0) {        /* rotate X    */
  323.     sine = -dy / hyp1;
  324.     cosine = hyp /hyp1;
  325.     } else {
  326.     sine = 0;
  327.     cosine = 1.0;
  328.     }
  329.     mat[1][1] = cosine;
  330.     mat[1][2] = sine;
  331.     mat[2][1] = -sine;
  332.     mat[2][2] = cosine;
  333.     i_multmatrix(mat);
  334.     mat[1][1] = mat[2][2] = 1.0;    /* be careful here to reinit    */
  335.     mat[1][2] =    mat[2][1] = 0.0;    /* those modified by the last    */
  336.                     /* paragraph    */
  337.     if (hyp != 0.0) {            /* rotate Y    */
  338.     sine = dx / hyp;
  339.     cosine = -dz / hyp;
  340.     } else {
  341.     sine = 0;
  342.     cosine = 1.0;
  343.     }
  344.     mat[0][0] = cosine;
  345.     mat[0][2] = -sine;
  346.     mat[2][0] = sine;
  347.     mat[2][2] = cosine;
  348.     i_multmatrix(mat);
  349.     i_translate(-vx,-vy,-vz);    /* translate viewpoint to origin */
  350. }
  351.  
  352. /*
  353.  *     window - 
  354.  *        Set up a 3d PERSPECTIVE window. This loads the matrix onto
  355.  *     the stack destroying the current transformation.
  356.  */
  357. i_window(left, right, bottom, top, near, far)
  358. float left, right, bottom, top, near, far;
  359. {
  360.     float Xdelta, Ydelta, Zdelta;
  361.     float mat[4][4];
  362.  
  363.     Xdelta = right - left;
  364.     Ydelta = top - bottom;
  365.     Zdelta = far - near;
  366.     if (Xdelta == 0.0 || Ydelta == 0.0 || Zdelta == 0.0) {
  367.     fprintf(stderr,"window: window width, height, or depth is 0.0\n");
  368.     return;
  369.     }
  370.     mat[0][0] = near * 2.0/Xdelta;
  371.     mat[1][1] = near * 2.0/Ydelta;
  372.     mat[2][0] = (right + left)/Xdelta;        /* note: negate Z    */
  373.     mat[2][1] = (top + bottom)/Ydelta;
  374.     mat[2][2] = -(far + near)/Zdelta;
  375.     mat[2][3] = -1.0;
  376.     mat[3][2] = -2.0 * near * far / Zdelta;
  377.     mat[0][1] = mat[0][2] = mat[0][3] =
  378.     mat[1][0] = mat[1][2] = mat[1][3] =
  379.     mat[3][0] = mat[3][1] = mat[3][3] = 0.0;
  380.     i_loadmatrix(mat);
  381. }
  382.  
  383. /*
  384.  *     perspective -
  385.  *        Set up a perspective window. This loads the matrix onto 
  386.  *    the stack destroying the current transformation.
  387.  */
  388. i_perspective(fovy,aspect,near,far)
  389. float fovy, aspect, near, far;
  390. {
  391.     float Zdelta, cotangent;
  392.     float matrix[4][4];
  393.  
  394.     if (fovy <= 0.1 || fovy > 180.0) {
  395.      fprintf(stderr,"perspective: fovy is out of range [0.1,180.0]\n");
  396.     return;
  397.     }
  398.     Zdelta = far-near;
  399.     if (aspect == 0.0 || Zdelta == 0.0) {
  400.     fprintf(stderr,"perspective: window width or depth is 0.0\n");
  401.     return;
  402.     }
  403.     fovy = ((fovy*M_PI)/180.0); 
  404.     fovy /= 2;                /* take half the angle*/
  405.     i_initmat(matrix);
  406.     cotangent = cos(fovy) / sin(fovy);
  407.     matrix[0][0] = cotangent / aspect;
  408.     matrix[1][1] = cotangent;
  409.     matrix[2][2] = -(far+near)/Zdelta;
  410.     matrix[2][3] = -1.0;
  411.     matrix[3][2] = -(2.0*far*near)/Zdelta;
  412.     matrix[3][3] = 0.0;
  413.     i_loadmatrix(matrix);
  414. }
  415.  
  416. /*
  417.  *     ortho -
  418.  *        Set up a 3d window. This loads the matrix onto the stack
  419.  *     destroying the current transformation.
  420.  */
  421. i_ortho(left, right, bottom, top, near, far)
  422. float left, right, bottom, top, near, far;
  423. {
  424.     float Xdelta, Ydelta, Zdelta;
  425.     float matrix[4][4];
  426.  
  427.     Xdelta = right - left;
  428.     Ydelta = top - bottom;
  429.     Zdelta = far - near;
  430.     if (Xdelta == 0.0 || Ydelta == 0.0 || Zdelta == 0.0) {
  431.     fprintf(stderr,"ortho: window width, height, or depth is 0.0\n");
  432.     return;
  433.     }
  434.     i_initmat(matrix);
  435.     matrix[0][0] = 2.0/Xdelta;
  436.     matrix[3][0] = -(right + left)/Xdelta;
  437.     matrix[1][1] = 2.0/Ydelta;
  438.     matrix[3][1] = -(top + bottom)/Ydelta;
  439.     matrix[2][2] = -2.0/Zdelta;        /* note: negate Z    */
  440.     matrix[3][2] = -(far + near)/Zdelta;
  441.     i_loadmatrix(matrix);
  442. }
  443.  
  444. /*
  445.  *     ortho2 -
  446.  *        Set up a 2d window. This loads the matrix onto the stack
  447.  *     destroying the current transformation.  Note that the new 
  448.  *    transformation leaves Z negated.
  449.  */
  450. i_ortho2(left, right, bottom, top)
  451. float left, right, bottom, top;
  452. {
  453.     float Xdelta, Ydelta;
  454.     float matrix[4][4];
  455.  
  456.     Xdelta = right - left;
  457.     Ydelta = top - bottom;
  458.     if (Xdelta == 0.0 || Ydelta == 0.0) {
  459.     fprintf(stderr,"ortho2: window width or height is 0.0\n");
  460.     return;
  461.     }
  462.     i_initmat(matrix);
  463.     matrix[0][0] = 2.0/Xdelta;
  464.     matrix[3][0] = -(right + left)/Xdelta;
  465.     matrix[1][1] = 2.0/Ydelta;
  466.     matrix[3][1] = -(top + bottom)/Ydelta;
  467.     matrix[2][2] = -1.0;
  468.     i_loadmatrix(matrix);
  469. }
  470.  
  471. float disttobezier(x, y, v1, v2, v3, v4)
  472. float x, y;
  473. float v1[2], v2[2], v3[2], v4[2];
  474. {
  475.     float vmin, vv;
  476.     float x1[2], x2[2], x3[2], x4[2], p[2], t;
  477.     float t3, t2, t1, t0;
  478.  
  479.     x1[0] = v1[0]; x1[1] = v1[1];
  480.     x2[0] = v2[0]; x2[1] = v2[1];
  481.     x3[0] = v3[0]; x3[1] = v3[1];
  482.     x4[0] = v4[0]; x4[1] = v4[1];
  483.     vmin = 1.0e30;
  484.     for (t = 0.0; t <= 1.0001; t += .001) {
  485.     t3 = t*t*t;
  486.     t2 = 3.0*t*t*(1.0-t);
  487.     t1 = 3.0*t*(1.0-t)*(1.0-t);
  488.     t0 = (1.0-t)*(1.0-t)*(1.0-t);
  489.     p[0] = t0*x1[0]+t1*x2[0]+t2*x3[0]+t3*x4[0];
  490.     p[1] = t0*x1[1]+t1*x2[1]+t2*x3[1]+t3*x4[1];
  491.     vv = (p[0]-x)*(p[0]-x) + (p[1]-y)*(p[1]-y);
  492.     if (vv < vmin) vmin = vv;
  493.     }
  494.     return sqrt(vmin);
  495. }
  496.